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The practical benefits of the hyper-accuracy properties of the discontinuous Galerkin 
method are examined. In particular, we demonstrate that some flow attributes exhibit 
super-convergence even in the absence of any post-procesing technique. Theoretical anal- 
ysis suggest that flow features that are dominated by global propogation speeds and decay 
or growth rates should be super-convergent. Several discrete forms of the discontinu- 
ous Galerkin method are applied to the simulation of unsteady viscous flow over a two- 
dimensional cylinder. Convergence of the period of the naturally occurring oscillation is 
examined and shown to converge at 2p+l, where p is the polynomial degree of the discon- 
tinuous Galerkin basis. Comparisons are made between the different discretizations and 
with theoretical analysis. 


I. Introduction 

The discontinuous Galerkin (DG) method provides a rigorous and robust means of formulating high-orcler 
methods for non-smooth unstructured grids. Such a method can be used to perform highly accurate simula- 
tions about complex geometries using grids that are readily generated. However, all high-order methods are 
more expensive than conventional second-order methods, on a per degree of freedom basis, so it is important 
to maximize the benefit of their high-order properties. The DG method is a projection technique in which the 
solution in each individual element is represented as an expansion in a local basis set; generally polynomials 
of degree < p. Formal proofs of error convergence 1-3 predict that the DG method with a basis of degree 
p should produce solutions that converge at a rate of p + a on general grids, and at p + 1 on non-smooth 
Cartesian grids and in several other special cases. Numerical experiments usually produce a convergence 
rate of p + 1 even on most non-smooth unstructured grids. 

However, there is a significant collection of analyses and numerical evidence indicating that the DG 
method provides a super-convergent solution of order 2p + 1 or higher. 4-10 In particular, a spatial eigenvalue 
analysis by Hu and Atkins 7, 8 for scalar convection and acoustic propagation clearly shows that the amplitude 
and phase errors converge respectively at rates of 2p + 1 and 2p + 2. The analysis demonstrated the super- 
convergence property for a range of p , and it was conjectured that it held for all p. This conjecture was later 
proven in Ref. 9. In Ref. 10, a spatial Fourier analysis of several variations of DG applied to diffusion for 
p = 4 predicted the convergence of the temporal eigenvalue to be p = 9 or p = 12 depending on the form of 
the discretization. It has been conjectured that perhaps different DG discretizations converged at rates of 
2p + 1 or 3 p. However, it will be shown in this work that this is not the case. 

Several post-processing methods have been developed and demonstrated for model problems. 6, 11-16 Gen- 
erally, these methods formally construct a solution of degree 2 p from a simulation of degree p by computing 
the convolution of the solution against a specialized kernel function. However, there are several drawbacks 
to these techniques that limit their usefulness. Although the kernel is local and formally compact, the kernel 
stencil size is proportional to p, and therefore it can be quite large. Yet, it has been shown that the stencil 
must be small relative to any wavelength of interest in order to improve that mode of the solution. 16 These 
methods usually require special grids that are difficult to apply to realistic flow simulations about complex 
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geometries. Finally, it is shown in Ref. 16 that some numerical errors, such as those produced by approxi- 
mate numerical boundary conditions, can introduce errors that the post-processing kernel cannot distinguish 
from the correct physical solution. Consequently, such numerical errors cannot be removed by any local or 
compact post-processing technique. 

Given that most post-processing techniques are difficult to successfully apply to a realistic simulation, 
it is reasonable to question whether the super-convergence property of the DG method is of any real value. 
As will be shown, the answer is yes; super-convergence is of real value. This paper examines the practical 
benefits of the super-convergence properties of the DG method. In particular, we demonstrate that certain 
flow attributes exhibit super-convergence in the absence of any post-procesing technique. Several discrete 
forms of the DG method are applied to the simulation of unsteady viscous flow over a two-dimensional 
cylinder. Convergence of the period of the naturally occurring oscillation is examined and compared with 
theoretical predictions. 

The first section describes the implementation of DG to the Navier-Stokes equations. The second section 
reviews the superconvergence properties of the DG method for advection and diffusion equations, and presents 
additional new analysis for diffusion equations. The third describes the numerical experiments performed 
and their results. 


II. Governing Equations and DG Formulation 

A. Navier-Stokes Equations 

The non-dimensional compressible Navier-Stokes equations in conservation form are: 


dp 

dt 

dpe 

~dt 

dpui 

~df 


9(pv-j) 

dxi 


= 0, 


d (huj — UjTjj + qj) 

dxj 

d (pujUj + Sjjp - Tjj) 

dxj 


= 0, 


= 0, 


(1) 

(2) 

(3) 


where p is the density, e is the internal energy per unit mass, iq is the component of velocity in the cartesian 
coordinate direction Xi, 

h = pe + p, 


Qi = 


— PrRe 

dxi 


Tij = pRe, 


-l 


7-1 


dui duj 2 duk 

dxj dxj 1 J 3 dxk 


T = p/p = (7 — l)(e — UkUk/ 2), Pr is the Prandtl number, and Re r is the Reynolds number based on 
reference state of the non-dimensionalization. The length and the thermodynamic variables have been non- 
dimensionalized with respect to a general reference state: 27 = Xi/L r , p = p/p r , p = p/p r , T = T/T r , and 
p = p/ p r where ' denotes climemsional quantities, and the subscript r denotes the reference state. It follows 
that u = u/u r , t = t/{L r /u r ), e = e r /u/, and Re r = p r u r L r /p r where u r = \Jp r / Pr- 

DG is applied to each of the equations 1-3 in essentially the same manner; however, the inviscid and 
viscous terms are treated differently. To facilitate the following discussion, each equation is cast in the 
general form: 

^ + V • (Fi - F„) = 0, (4) 

where U denotes either p, pe, or pu.i , and F, and F„ denote corresponding the inviscid and viscous contribu- 
tions to the flux. Equation (4) is transformed to a local computational coordinate system for each discrete 
element. Let (£,77, () denote the local coordinate system with Jacobian J = d(x±, X 2 , x^j/d^, rj, Q. In the 
local coordinates, equation (4) has the same form: 


dU 

— + V • (F, ; - F„) = 0, 


( 5 ) 


but with U = U\J\, Ft = J|, and F v = J^t^J |. 
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B. DG Discretization 

The computational domain is subdivided into non-overlapping elements that cover the domain. The DG 
discretization is formulated locally in each element in a similar manner. The solution within each element is 
approximated as an expansion in a local basis set {&&}, usually polynomials of degree < p, 

N 

U ^ ' b k U k ^ 

k=0 


where N denotes the number of terms in the basis set. The current work will uses two-dimensional monomials 
of the form f for all (i,j) pairs satisfying 0 < i + j < p. Following the quadrature-free formulation, 17 the 
fluxes are expanded in a similar manner, 

M M 

F^ — ^ ' bA, k and F„ — ^ ) b k ^v,kt 
k — 0 k = 0 

however, the degree of the flux expansions are allowed to be higher than that of the solution. The number 
of unknowns in each element is the number of physical variables times the size of the basis set, N. An equal 
number of equations governing these unknowns is derived by multiplying the governing equations by each 
member of the basis set, and integrating over the element. The integrals of the flux terms are integrated by 
parts to obtain the following weak form: 

[ bk~^~d£l — [ V6fe(Fj — F„) dCl + [ b k (Fi - F„) • n ds = 0 Vfc : 0 < k < N, 

Jn ot Jq Jqq 

where fi denotes the element, dfl denotes the element boundary (or edge) and n denotes the outward unit 
normal. Because neighboring elements have independent local approximations for the solution and fluxes, 
the flux terms in the edge integrals are multi-valued. These fluxes are replaced with numerical fluxes that 
are functions of both the local and neighboring solutions. 

N N 

= Fi — ^ ( bkfi )k and F.y • = F v — ^ ) b k f v k . 

k = o k = o 


Because all of the fluxes are represented as expansions in the basis set, the integrations can be performed 
directly to produce a matrix equation of the form: 


M [wfe] + V • [f^fc — f^fc] + 

edges 


B, 


f *. k ~ 


= 0 


where M, V and B are derived in Ref. 17. 

Gradients of the solution required for the viscous terms are evaluated using the DG methodology in a 
similar manner. Let 

N 


<J = 


y b k <Tk = Vro, 
fc= o 


where w denotes either u,; or T. It immedately follows that 


b k crdn — / \7b k wdn+ / b k wnds = 0. 


I an 


As before, the multi-valued edge flux is replaced by a numerical flux, wn| an = um. 

The inviscid edge flux is modeled by an approximate Riemann flux. The current work uses a Lax-Fredricks 
flux of the form: 

F i = {F i }-n-\[U]/2 

where {x} denotes the average of x, [a:] denotes the jump in x (neighboring value minus local value), and A is 
the maximum absolute eigenvalue of the Jacobian of the system of inviscid fluxes. Quantities not explicitly 
inside an average or jump operator ({x} or [a:]) are evaluated using data from the local element. Note: 
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because the outward normal, n, is of the opposite sign in neighboring elements on a common edge, we have 
the identity: {F.;} • n = — [F* • n]. 

Three variations of DG for the discretization of the Navier-Stokes equations are evaluated in this work. 
They differ only in the treatment of the numerical viscous fluxes as defined in Table 1. The interior- 
penalty formulation, 20 denoted DG-IP, and the closely related method of Bassi 21 that uses the “lifting” 
operator to evaluate a, are the most commonly used forms of DG for Navier-Stokes simulations. The DG-IP 
discretization is compact, and it is perhaps the easiest in which to formulate exact Jacobians for use in 
implicit solvers and adjoint methods. The second and third forms are variations of the local-discontinuous 
Galerkin (LDG) method of Cockburn and Shu. 18 The second method, obtained by setting (3 = 0 (in Table 
1), uses a symmetric edge flux in both the gradient and physical flux terms, and also adds a penalty term to 
the physical flux. This method will be denoted LDG-C (for LDG-central). The LDG-C spatial discretization 
is also easy to implement, but it results in a non-compact stencil which complicates the implementation of 
both implicit solvers and adjoint methods. The final method uses “one-sided” operators for the gradient and 
physical flux terms that alternate direction between the two. 18,19 Hence, it is denoted as LDG-OS. It also 
results in a compact discretization; however, its implementation is slightly more complex than that of the 
DG-IP method. It is obtained by setting (3 = ±1/2 and a = 0. The sign of [3 is determined by imposing 
a global, but otherwise nearly arbitrary, ordering on the elements, and then choosing the sign according to 
whether the local element is “above” or “below” its neighbor. 


C. Non-linear Flux Terms 

The quadrature-free form of DG requires that all fluxes be expanded in the basis set. This process is trivial 
and exact when the equations are linear, but must be approximated for most non-linear cases. Previously, 22,23 
it was noted that the inviscid fluxes of the Euler equations consist of ratios of polynomials such as 

j = jpui){puj) or j = jpui){pe) 

P ' P 

Polynomial approximations to the flux can be forumulated by multiplying by p, and projecting back into 
the polynomial space. For example: 

J hpf = J b k (pUi)(pUj). 

Considering only design order arguments, the polynomial products can be truncated to degree p, and the 
projection can be limited to k < N. However, in Ref. 23 it was shown that extending the degree of the flux 
improved robustness in cases of shock capturing. In the current work, polynomial products are truncated to 
p+ 1, and the projection is over that same set of functions. The operation is efficiently performed by noting 
that the left-hand side is the same for all terms, 


LHS = 



allowing the LU - decomposition of the matrix to be re-used many times. 

The Uj and T terms, needed to compute the viscous gradients, are obtained in a similar manner by 
solving: 


J b k pf = R , where R = 


b k (pui) or 


/ 


b k P. 


Scheme 

w 

Tv 

LDG-C and LDG-OS 18,19 
DG-IP 20 

M - /3 • M 
M 

{fv} + /?[/„] - a Ml 
{Vw} - a(H) 


Table 1. Numerical viscous fluxes. 
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III. Super-convergence properties of DG 


For a DG solution that is an expansion in a finite polynomial basis, say degree = p, any direct discrete 
comparison with an exact solution will clearly result in a measured error that converges no faster than p+ 1 . 
It would seem that some post-processing would be necessary in order to obtain any real benefit. However, 
previous spatial eigenvalue/eigenfunction analysis of acoustic propagation ' -8 have definitively demonstrated 
that the eigenvalues of the DG solution are super-convergent. In its simplest form, this analysis assumes 
that the discrete eigensolution is of the form 

<(£, t) = e-^XMO (6) 

where A = e lKh and ■u(^) is a polynomial function of degree p within a single element at n = 0. Substitution 
into a model wave equation, with wave speed a, produces a dispersion relation that governs the eigensolution. 
For the compact DG discretization, this equation has the form: 

(1 - 7 )[G(LA)A - H(iK)\ + (-l) p+1 (l + _ H(-iK)] = 0 (7) 

A 

where K = — , and H and G are polynomial functions of the normalized frequency, K, and 7 arises out of a 
generalized representation of several flux splitting formulations, as described in Ref. 7. A conventional Fourier 
analysis sets the wave number, AA, and solves for the resulting temporal behavior. However, this results in 
an eigenvalue problem of size p + 1 that must be solved numerically, and that produces multiple solutions 
that require some interpretation. As noted in Ref. 7, equation (7) gives a simple quadratic relation for A 
that can be solved analytically when u> is given. One root corresponds to a non-physical mode that has been 
shown to be damped such that it does not interfere with the properties of the asymptotic physical solution. 
It was also shown that the ratio of polynomials, H(iK)/G(iK), equals the exact Pade approximation of e lK 
to order 2 p + 2 ; and consequently, the physically significant root is of order: 

A = e iK + O((coh) 2p+2 . 

Thus, the dispersion and dissipation errors converge as Re(iKh ) ~ K + O((ojh) 2p+3 ) and Im(Kh) ~ 
O((coh) 2p+2 ), respectively. Note that because the analysis is for a single time step (or cell translation), 
the convergence rate of the error after integrating to a fixed time is one order lower. 

In Ref. 10, a Fourier analysis of DG applied to diffusion indicated that all three discretizations described 
in the previous section have a super-convergence property. In that work, the analysis was performed for 
p = 4 only, and it was found that the DG-IP discretization converged at a rate of p = 9, while the LDG-C 
and LDG-OS discretizations converged at a rate of p = 12. It has been conjectured that this implied general 
rates of 2p + 1 for DG-IP and 3 p for LDG-C and LDG-OS. 

Here the earlier analysis is extended to include p = 1,2, and 3 in order to determine the true convergence 
properties. The procedure is similar to the analysis described above (ref. 7) except that here we choose Kh 
and solve governing equations for K. After including the diffusion terms, the governing equations define 
a matrix eigenvalue problem of size p+ 1 that must be solved numerically for the eigenvalues, AT*,. Each 
of the multiple solutions for obtained at a given Kf, correspond to spatial wave number that is shifted 
by some multiple of tt (e.g. Kh + kn). Rigorous association of Kj. with the appropriate wave number 
would require examining the wavelength of the eigenfunction. However, for the current purpose of inferring 
order properties, it is generally sufficient to simply sort the solutions from smallest to largest, (A) < Kj 
if i < j ) and associate each solution with a wave number as [AT*, Kh + kn}. Figures 1 a-c show the result 
of this analysis for each of the three discretizations, and for p = 1 — 4. The dashed line provides reference 
convergence rates. The curves are labeled p = a(b) where a is the degree of the basis set, and b is the 
convergence rate of the reference line. From these results it is clear that DG-IP and LDG-C converges at 
2p+2, and LDG-OS converges at 2p+4. As before, these rates are for a single time step, and the convergence 
rates at a fixed time would be one order lower (assuming St oc h). 

Both analyses predict a super-convergent relationship between the temporal and spatial eigenmodes (i.e. 
between e~ lut and A), but say nothing about the function u{£f). In fact, u{f ) is simply a polynomial function 
of degree p that is completely determined by initial and boundary conditions; and than evolves in time 
and space as governed by e~ lut X n . Post-processing techniques are generally local and compact, and can be 
used to construct a higher-order local solution. However, because post-processors are local and compact, 
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(a) DG-IP (b) LDG-C (c) LDG-OS 


Figure 1. Convergence rate predicted by Fourier Analysis 


they cannot correct the global dispersion and dissipation errors. The global super-convergence properties are 
inherent to the DG discretization method, and are a prerequisite for any successful post-processing technique. 
In essence, local post-processing techniques can only improve w(£) by “smoothing” high-order noise in the 
local solution. 

This analysis has several implications for realistic simulations. First, it implies that the accumulation 
of error remains small for features that evolve for long periods of time or space, perhaps smaller than the 
initial and boundary error. Thus for the scatter of acoustic waves off of a complex geometry, one could 
expect super-convergence of the wave propagation and the resulting interference pattern. Also, the growth 
of transitional and weak instability modes that evolve for long times and distances may be super-convergent. 

Second, this analysis suggests that naturally occurring phenomena that depend weakly on initial and 
boundary conditions may be super-convergent. For example, the profile of a boundary layer on a flat 
plate may be super-convergent; especially since the shape of the geometry and the constant pressure outer 
condition can be represented exactly by even the lowest order DG basis. Viscous layers over more complex, 
yet sufficiently smooth, geometries may also be super-convergent. Finally, naturally occurring unsteadiness 
such as that in bluff body and cavity flows, as well as transitional instabilities, may be super-convergent. 

IV. Test Case 


A. Model Problem 

Viscous flow over a two-dimensional cylinder is chosen as a test case because the flow produces a well behaved 
natural unsteadiness at a dominant frequency. In addition, the geometry is simple and analytic, but non- 
trival (e.g. it cannot be represented exactly by the polynomial basis.) Previous investigations concerning 
the super-convergence of acoustic scatter 10 suggest this boundary shape is sufficiently smooth so as not to 
obscure any super-convergence behavior that may exist. 

Figures 2a-c illustrate typical grids showing the full domain on the coarsest grid, and the coarsest and 
finest grids near the cylinder. Simulations were performed on a sequence of four grids with the number of 
elements around the half cylinder, Nx, equal to 5, 10, 20, and 40, for p < 3; and with Nx = 9, 12, 18, and 27 
for p = 4. Although the grids shown have straight edges, element edges lying on the cylinder are curved and 
defined by projections onto a polynomial space. Points 1, 2, and 3, shown in figure 2(b), are the locations 
at which the flow properties were sampled. 

The freestream Mach number, is 0.2 and the Reynolds number, Re , is 200, based on the cylinder 
diameter, d , and freestream velocity. Taking the reference state as freestream, the reference Reynolds number 
is given by, 
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= Re^j/i^Moo). 

The non-dimensional viscosity, /z, is held constant at 1. 

The simulations are started impulsively from a modified freestream. Near the cylinder, the flow velocity 
is reduced linearly in the radial direction such that the initial solution satisfies the no-slip condition. Addi- 
tionally, a compact, low amplitude vortical component is added to the initial solution in the neighborhood 
of the cylinder to excite the natural oscillation, and reduce the cost of the simulation. 


B. Error and Convergence Measurement Techniques 

The truncation error of any simulation is usually assumed to behave as 

OO 

Uk — U ex act = ot n {h k / ho ) n , (8) 

n=r 


where k identifies a grid within a consistent family of grids. For sufficiently small h k /ho, the asymptotic rate 
of convergence is given by 


n 


r k = log 


Uk U exact 

Uk-i — U ex act. 


/logOfc/ft.fc— i). 


(9) 


Of course in most realistic cases, H exac t is not known. It is common to use a very fine grid solution as 
a surrogate for the exact solution; however, this can be prohibitively expensive. A practical alternative is 
simply to note that the expression 


OO 

U k - U k+l = Y, Ml - (h k+1 /h k ) n )(h k /h 0 ) n , (10) 

n=r 


converges at the same asymptotic rate as true error given in equation (8). Thus, another approximation for 
the asymptotic rate is given by 


r k = log 


U k — U k+ 1 
U k - 1 — U k _ 


/\og{h k /h k -i). 


(11) 


The common practice of simply “fitting” the asymptotic error equation to the solutions on a sequence of 
three grids results in the same formula. This approximation for the convergence rate is accurate only if 
all three grids used in the evaluation of equation (11) are within the asymptotic regime. Note that non- 
monotonic convergence, (U k — U k+ \)/(U k -i ~ U k ) < 0, is a clear indicator that the coarser grids are not in 
the asymptotic regime. 

A third approach, described in Ref. 24 and 25, uses a sequence of 4 grids to relax this requirement by 
modeling the error as 

A*, = U k — U k+ \ = a r (h k /ho) r + a r+ i(h k /ho)' +1 . (12) 

The convergence rate is approximated by 


r k = log 2 


3Afc + i If 3Afc + i \ 
4Afc + 2 V \4Afc +2 / 


A k 

2Afc + 2 


where h k -\/h k = 2 Vfc. (13) 


Here, this formula is generalized to arbitrary mesh ratios, a = h k -\/h k as: 


r k = log,j 


rAfc + i //TAfc + i\" Afc 
Afc+2 V \ Afc + 2 ) aA k+ 2 


(14) 


where T = (a + 1) / (2 a). If all 4 grids are in or near the asymptotic regime, then A k+ i/ A k+2 > 0, and “+” 
is the appropriate sign before the radical. References 24 and 25 give a quantitate means of checking the 
validity of equation (14) by computing the ratio of the high- and low-order terms in equation (12). 


a rhk = 4(1 — 2 P ) (A fc — 2 p A fc+1 ) 
a r+1 (1 - 2 P+ 1 ) (2P+!A fc+1 - A fc ) ' 
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If this term is small (< 0(1)), then it is likely that all 4 grids are within or near the asymptotic regime. 
However, if it is much greater than 1, then it is likely that higher-order terms neglected by equation (12) 
are contributing to the error on the coarser grids, and that equation (14) is unreliable. The formulas given 
in equations 11 and 14 will be referred to in the remaining discussion as 3-grid and 4-gricl approximations, 
respectively. 


V. Numerical Test 

The three variations of DG described above are applied to simulate unsteady viscous flow over a two- 
dimensional cylinder. For each discretization method, the solution is approximated by a polynomial basis of 
degree p = 1,2,3, or 4, which should result in a design order of accuracy of p + 1 = 2,3,4, or 5, respectively. 
Simulations on a sequence of 4 grids were performed for each discretization method and for each basis degree. 

All flow solution variables were monitored at the three sample points shown in figure 2(b). The value of 
pv at point 3 will be used to illustrate typical results in the remaining discussion. Figure 3 compares the 
evolution of the p = 1 solution on the coarsest grid with the p = 3 solution on the finest grid. Both solutions 
experience a brief transient, but settle to a quasi-periodic solution after 6-8 cycles. There is a significant 
difference in the period between the two solutions: approximately 30%. 



Figure 3. Evolution of pv component at point 3. 


The period, A, is computed as the time interval between events when the solution crosses a specified 
value. The specified value varies with each sample point location and flow variable. Even though the time 
step is extremely small due to the explicit time marching method used, a high order interpolation is used 
to determine the precise moment of this crossing event between time steps. Before comparisons can be 
made between solutions on different grids, it is necessary confirm that transients in the solutions are small 
relative to the differences measured between grids. Figure 4 provides a measure of the decay of non-periodic 
transients by showing the change in the period between successive periods for all grids, basis degrees, and 
discretizations. Results for the different discretizations are shown without distinction, and generally cannot 
be distinguished from one another except at isolated points. Generally speaking, transients in all grids, 
degrees, and discretizations decay at similar rates, but the more accurate solutions are less monotone in 
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their convergence. 



Figure 4. Decay of non-periodic transients. Results from different discretizations are shown without distinc- 
tion. 

The change in the period with mesh refinement, A*, — Afc_i, is shown in figure 5 for all grids, degrees, 
and discretizations. Again, the different discretizations are shown without distinction. The higher-order 
cases require 9-11 periods before this metric begins to reach stable values. This is expected because the flow 
transients, shown in figure 4, must decay to levels well below the change in period with mesh refinement. 
After the initial transients decay, the results from the different discretizations overlay one another. Later in 
this discussion, subsets of this collection illustrating typical results for single discretization will be isolated 
and examined separately. 

Grid convergence of A k — Afc-i is shown in figures 5. Because all three discretizations give similar results 
for p < 3, only for LDG-OS computations were performed at p = 4. Recall that A& — A&_i converges 
asymptotically at the same rate as the true error. Also note that if the error on a sequence of grids is in the 
asymptotic regime, then the curves will be uniformly spaced when viewed on a log plot. However, only the 
p = 1 cases (figure 5a) display this behavior. This indicates that all four grids are within the asymptotic 
regime for all cases with degree p = 1. Also, note that for p = 3 (figure 5c) the change in the period is 
small for the coarse grid cases ( Nx = 10). This does not mean that the error in the period is small; only 
that the change in the period between the two coarsest grids is small. This generally indicates that the 
convergence is not monotone, and that the coarsest grids are well outside of the asymptotic regime for p = 3. 
Directly examining the data for p = 2 (figure 5b) reveals that the convergence is highly oscillatory. The sign 
of the change in the period flips between every mesh refinement (positive for Nx = 10 and 40, negative for 
Nx = 20). This suggest that all of the grids may be well outside of the asymptotic regime, and that more 
than two error modes are contributing to the error on most, if not all, grids. Therefore, neither the 3-grid 
or the 4-grid approximation for the convergence rate can be applied to the p = 2 case. 

Convergence rates for LDG-OS, computed using both 3-grid and 4-grid approximations, are shown in 
figure 6. Cases with p = 1,3, and 4 super-converge at a rate approaching 2p + 1 that is steady beyond 
the 11th cycle ( t > 400). The 3-grid approximation tends to be lower than the 4-grid approximation by an 
amount that increases with p. This is probably because the higher order discretizations require finer grids to 
be firmly in the asymptotic regime, and the 3-grid approximation is more strongly influenced by neglected 
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(a) p = 1 


(b) P = 2 




(c) p = 3 


(d) p = 4 


Figure 5. Change in period with mesh refinement. 
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terms. Because the p = 1 case is converging uniformly on all 4 grids, the 3-grid approximation can be applied 
to both the three coarser grids, and the three finer grids. The two computations give the similar results 
and both results are shown without distinction (the solid red squares.) Again, neither approximation for the 
rate of convergence is directly applicable to the p = 2 case; however, the absolute levels of the change in the 
period are similar to those of the p = 3 case. 

Although there are theoretical reasons to expect that LDG-OS might exhibit a super-convergent rate 
of 2 p + 3, it was not observed in these numerical experiments. This may be because the inviscid terms, 
which are treated the same in all cases, are the limiting contribution to the error. It is known that the 
Strouhal number is essentially constant over a range of Reynolds and Mach numbers, which means the 
shedding frequency is essentially proportional to the flow velocity. Although the viscous layer plays a vital 
roll in vortex formation, it is likely that the frequency of shedding is strongly influenced by the rate at which 
vortices advect downstream, and acoustic feedback travels upstream to the cylinder. 
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Figure 6. Convergence rate of LDG-C predicted by 3-grid and 4-grid approximations. 


VI. Conclusions 

The several discrete forms of the discontinuous Galerkin method have been applied to an unsteady viscous 
flow to evaluate the convergence rates of naturally occurring oscillations. Numerical tests confirm that the 
super-convergence predicted by analysis can be realized for all three discretizations. The period computed by 
all discrete forms of discontinuous Galerkin for p = 1,3 and p = 4 converge at a rate approaching 2p+ 1. The 
change in period with mesh refinement oscillates in sign for the p = 2 case, and formulas for the convergence 
rate do not directly apply. Although the change in period with mesh refinement for p = 2 is similar to that 
of p = 3, the error of all p = 2 cases is probably much larger than that of p = 3. These results clearly 
demonstrate that important flow features from simulations using the discontinuous Galerkin method are 
super-convergent without the application of any post-processing technique. 
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